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Abstract 

A numerical method is developed for calculating the 
real time Green's functions of very large sparse Hamil- 
tonian matrices, which exploits the numerical solu- 
tion of the inhomogeneous time-dependent Schrodinger 
equation. The method has a clear-cut structure reflect- 
ing the most naive definition of the Green's functions, 
and is very suitable to parallel and vector supercom- 
puters. The effectiveness of the method is illustrated 
by applying it to simple lattice models. 

1. Introduction 

In many fields of quantum physics, evaluation of 
the Green's functions constitutes the most important 
and difficult part of the theoretical treatment"'^-'. For 
example, to compute physical quantities of manybody 
systems in condensed matter physics, one should often 
calculate the Green's function of Hamiltonian matrices 
having a degree N of 10^ or more. Therefore efficient 
numerical algorithms, such as recursive Green's func- 
tion methods, quantum Monte Carlo methods, and the 
Lanczos methods have been developed and applied to 
various problems. 

Recursive Green's function methods^^ have succeeded 
in evaluating dynamic quantities of relatively small 
systems by calculating directly the real-time Green's 
function. For example, the conductance of quantum 
dots in chaotic, and regular regimes has been inten- 
sively investigated with these methods^-'. However, 
this scheme is prohibitive for huge Hamiltonian matri- 
ces because it requires computational time increasing 
rapidly as a function of the matrix size. 

Quantum Monte Carlo methods^'^\ which gener- 
ate the imaginary-time Green's functions, have been 
successfully used for evaluating thermodynamic quan- 
tities of relatively large systems. For evaluating dy- 
namic quantities such as conductivity, however, one 
has to rely on numerical analytic continuation (e.g., 
maximum entropy methoc^^) from the imaginary-time 
Green's functions to the real-time ones. This proce- 
dure is, however, not unambiguous due to two reasons: 
one is the statistical errors originating from Monte 
Carlo sampling, which are amplified by numerical an- 
alytical continuation, and the other is the bias intro- 
duced by the default model in the maximum entropy 

t Condensed from the article submitted to Phys. Rev. E. 



method. 

The Lanczos methods^'®^ have been one of few re- 
liable techniques for evaluating dynamical responses of 
moderate-size Hamiltonian matrices. The Lanczos meth- 
ods use a linear transformation to a new basis in which 
the Hamiltonian matrix has a tridiagonal form, and 
lead to a continued fraction representation of the di- 
agonal matrix elements of the Green's function. The 
drawback of these methods is the numerical instability 
which may lead to spurious eigenstates^K Recently, 
the Lanczos method has been extended to the finite 
temperature case by introducing random sampling over 
the ground and excited states^"^. 

In this paper, we present a new algorithm called 
Particle Source Method (PSM), which is based on the 
most naive and effective definition of the real-time Green's 
functions^^^ . Namely, we calculate numerically the 
time-dependent Schrodinger equation having a source 
term, and see how the wave function responds to the 
particle source. 

This method has resemblances to the Forced Os- 
cillator Method (FOM), which has been developed by 
Williams and Maris and applied to various classical 
and quantum problems^'^'^'*-'. The FOM calculates the 
classical equations of motion of the coupled harmonic 
oscillators driven by a periodic external force, where 
the matrix elements of the Hermitian matrix give the 
frequency and the coupling of the fictitious oscillators. 
Our method is, however, much more clear-cut than 
their method, when applied to quantum systems, since 
we calculate the time-dependent Schrodinger equation 
itself instead of the classical equations of motion mapped 
from the quantum Hamiltonian matrix. The differ- 
ence between the two methods is analogous to the dif- 
ference between the old quantum theories describing 
electronic states of an atom as an ensemble of ficti- 
tious harmonic oscillators and the modern quantum 
mechanics describing them by bra's and ket's. 

Preceding to the present article, several authors 
have already solved the homogeneous and inhomoge- 
neous time-dependent Schrodinger equations numeri- 
cally I5'i6>i7,i8,i9) ^ Most of them are, however, inter- 
ested in launching wave packets in the computer and 
watching them move around. Several of them tried to 
extract time-independent quantities from the motion 
of the wave packets'^*'^^-'. Unfortunately, their inter- 
est was limited in obtaining several eigenvalues and 
eigenvectors of the Hamiltonian, but not the Green's 



functions. As the result, they could obtain only the 
exact peak position of the spectrum function (i.e., the 
imaginary part of the diagonal elements of the Green's 
function), but could not calculate the correct shape 
of the spectrum function, the real part, and the off- 
diagonal elements. This is in contrast to our method, 
which can calculate both real and imaginary parts, and 
both diagonal and off-diagonal elements of the Green's 
function without calculating the eigenvalues and eigen- 
vectors. 

In section we present the basic ideas of the PSM. 
In section we extend the PSM to the finite temper- 
ature case. In Section we present numerical exam- 
ples to illustrate the effectiveness of the methods. A 
summary is given in section [s]. 

2. Particle Source Method 

2.1. Single frequency calculation 

Let us introduce the time-dependent Schrodinger 
equation with a time-dependent source term. 



H\(t>-t) + |j)e-*('^+*'')*6»(i) 



(1) 



where the wave function \(j)]t) and an arbitrary source 
|j) are A^-component complex vectors, the Hamilto- 
nian H is an N x N Hermitian matrix, to is the fre- 
quency of the source, and 77 is a small positive imagi- 
nary part of the frequency, which determines the reso- 
lution of frequency. Note that this source term grows 
up exponentially as a function of time due to this small 
positive number, which simulates adiabatic switching 
on of the particle source. This adiabatic switching on, 
which has been absent in the preceding works^^''^'^'^^-', 
is essential to calculate the exact shape of the Green's 
function as a function of energy. 

The solution of this equation with the initial con- 
dition = 0) = becomes^^ 
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where we have neglected the second term in the paren- 
theses of (H). This approximation is justified by using 
sufficiently long time ti satisfying the condition 



(6) 



where S is the required relative accuracy of the Green's 
function. 



Then, from the Fourier transformation of (||), the 
Green's function operated on the ket \j) is obtained as 

^ dt' |0;t)e*(^+*'')* 



(7) 



= T- / dt'G{co + ^^)\j) 
Jo 



If only one or few matrix elements are necessary, we 
can calculate only these matrix elements as 

- / 'di'(i|0;t)e*('^+*'')*' 
ti Jo 

= 1 dt'{t\Gicj + iv)\j) 
f-i Jo 

= (^\G{w + i7^)\3)^G,,{uj + iii) (8) 

where (i| is an arbitrary bra. 

Since the numerical error due to the finite timestep 
is proportional to (wAt)'^ ^''^ the best choice of lo is 
a; = 0. The matrix elements with energy w ^ can be 
obtained by calculating the shifted Green's function at 
cj = 



G"(w = 0;77) = 



1 



O + i-q- H' 
with the shifted Hamiltonian 

H' = H -Lul 
where / is the unit matrix. 



(9) 



(10) 



2.2. Multiple frequency calculation 

Let us introduce the time-dependent Schrodinger 
equation with a multiple frequency source term. 



dt 



emu) 



where loi = lAuj. 

The solution of this equation with the initial con- 
dition \(f);t ^ Q) =0 becomes 
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where we have neglected the second term in the paren- 
theses of (|l^) as in the single frequency calculation. 



Then, from the Fourier transformation of ( |15D , the ma- 
trix elements of the Green's function are obtained as 

- / 10; t)e'('^'' 

h Jo 



'■2 Jo 



(i\G{uji> + if])\j) 
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where we have neglected the second term in . This 
approximation is justified by using sufficiently long 
time ^2 satisfying the condition 



t2Aw > 1/5. 



(19) 



2.3. Analysis of the Numerical Method 
2.3.1. Solving the Schrodinger Equation 

To solve the time-dependent Schrodinger equation 

(^ numerically, we discretize it by using the leap frog 
methodi5.i6,i7) 



;t + At) 
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-2iAt\j)e-'^'^+''^'>'e{t). (20) 



where At is the time step. The time step is set as 

At = a/E^ax (21) 

where Emax is the absolute value of the extreme eigen- 
value. We usually use the parameter a between 10~^ 
and 10"^ 

Another method for the time-dependent Schrodinger 
equation is the Suzuki- Trotter decomposition of the 
time-evolution operator. Though the Suzuki- Trotter 
decomposition can be applied effectively only to a spe- 
cial class of Hamiltonian, it might have the advan- 
tage of the leap frog method. First it allows larger 
time step. Second it can be used with non-Hermitian 
Hamiltonian, such as the Hamiltonian with absorbing 
boundary condition. 

2.3.2. CPU time 

The computational time to calculate G{LO + ir])\j) is 
estimated by the number Np^od of matrix-vector prod- 
ucts in (pO|), which is equal to the integration time t 
devided by time step At, 
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t 

At 
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(22) 



Introducing (||) into (p^), we obtain the number of 
matrix-vector products for the single frequency calcu- 
lation 



jyprod _ 



log 6 E„ 



a 



(23) 



Therefore the relative error S becomes exponentially 
small as a function of computational effort Nprod- On 
the other hand, the resolution r/ is inversely propor- 
tional to Nprod, that is, we need longer CPU time for 
higher resolution. 

Introducing (^ into (|2^), we obtain the number 
of matrix-vector products for the multiple frequency 
calculation 



od 



1 E„ 
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Therefore the relative error 6 is inversely proportional 
to Nprod, which means slower convergence than the sin- 
gle frequency calculation. The distance between the 
frequencies to be measured. Aw, is inversely propor- 
tional to Nprod, that is, we need longer CPU time as 
we increase the number of the frequencies. Actually, 
the integration time for the multiple frequency calcu- 
lation should be the longer than ti and t2. However, 
because ti < t2 is usually satisfied, t2 determines the 
CPU time for the multiple frequency calculation. Since 
the computational effort for a product of sparse ma- 
trix and vector is proportional to the matrix size N , 
the total computational time is estimated as 

CPU _ -^ogSE„ 



-Nm 



nCPU 



1 E„ 



V 
-N 



(25) 
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where N'^ is the number of the frequencies to be mea- 
sured. 

2.3.3. Calculation of G{uj - ii]) 

So far, we have been calculating the Green's func- 
tion whose frequency has a positive imaginary part. 
When we need to calculate G{uj — irj), we substitute t 
by — t, and rj by —rj in (|l|) an d fol low the same pro- 
cedure as described in section |2.1.| . Then we obtain 
G{uj - ir])\j). 

2.3.4- Product of the Green's functions 

Since \j) in (|l|) is an arbitrary ket, we can repeat 
the calculation of the Green's function by using a new 
source term, 

|j2)e-'('^^+*''^)*6'(i) = 

AM^i + *m)|j)e-'("^+"'^)*e(t) (27) 

where Ai is an arbitrary operator whose matrix el- 
ements are known. In general, we can calculate the 
matrix elements of a product involving several Green's 
functions and other operators as 

{i\AnG{ujn±iVn) ■ ■ ■ A2G{uj2±im)AiG{uji±irii)Aa\j).{28) 

2.3.5. Remote Eigenvalue Problem 

The remote eigenvalue problem pointed out in ref- 
erence 19) does not appear in our methods, since we 
use very small time step in order to integrate the Schrodinger 
equation stably by using the leap frog method. 



2.4- Application to Manybody Problems 
2.4-1- Single-particle Green's function 

We can apply our methods for calculating the Green's 
functions of an N-particle system at the ground state. 
As an example, let us see how we can calculate the re- 
tarded single-particle Green's function of an N electron 
system on a finite lattice, 

Gij (uj + irj) 

/oo 
dr(5|{a.(r),a](O)}|5)e'("+''')^0(r) 
-OO 

= +{g\a,G{Eg + uj + iri)a]\g) 



{g\a]G{Eg - u - ir])a,\g) 



(29) 



where and Oj are the annihilation operator at site i 
and the creation operator at site j; \g) and Eg are the 
groundstate of the N electron system and its energy, 
respectively. Since each term of ( p9| ) has the form of 
(p8|), we can calculate Q{uj) as follows: 

First, we calculate the ground state \g)N of the 
N-electron system by using one of existing methods 
such as the Lanczos method, the quantum Mote Carlo 
method, or the finite difference method^^^. The ad- 
vantage of using the finite difference method is that it 
can recycle most of the subroutine resource written for 
calculating the Green's function since both programs 
solve the time-dependent Schrodinger equation in the 
same way. 

Second, we operate to the ground state to ob- 
tain an A'^ -|- 1 electron state, |j)Ar+i = a]|ff)7V- In a 
similar way, we calculate another -I- 1 electron state, 

\i)N+i = al\g)N- 

Finally, the retarded Green's function is calculated 
in A'^-l-l-electron subspace using the method in the pre- 
vious subsection together with the state vectors | j)jv+i 
and \i)N+i- 

2.4.2. Optical Conductance 

The optical conductivity is expressed within the lin- 
ear response theory as 

a,j,^{LO -t- irj) 

= - / die^^-+^'')*(g|j.(t)j,(0)|.g) 

_2 

= — l^g\2^G{uJ + Eg+l^)j^\g), (30) 
which we can calculate by using the PSM. 

3. Monte Carlo Particle Source Method 

In this section, we extend PSM to finite tempera- 
ture case by using Monte Carlo Particle Source Method 
(MCPSM), a combination of PSM and the Monte Carlo 
method for calculating the trace of a large matrix. It 



turns out that, for sufficiently large systems, only a sin- 
gle configuration of random variables suffices to eval- 
uate the desired expectation value at finite tempera- 
tures. 

3.1. Monte Carlo Calculation of Trace 

Computing trace of a large matrix A requires evalu- 
ation of N diagonal elements of the matrix. Therefore 
it would take formidable computational time if we try 
to use PSM for calculating the trace of the product of 
operators including the Green's functions. A Monte 
Carlo method to estimate trace of large matrices^^-* 
makes it possible to evaluate the trace of this kind. 

Let us introduce a set of random variables , {n = 
1, • • • , N) that satisfy the relation 



)) = Sn'n 



(31) 



where (( • )) implies statistical average. Then we de- 
fine a random ket as 



N 



i'f) = E 



(32) 



with the chosen basis set {|n)}. Then the statistical 
average of ($| gives an approximation of the trace, 

(( » 



iv[A]. 



{n'\A\n) (33) 
(34) 



The second term in ( 33 ) gives the statistical error when 
the average is evaluated by the Monte Carlo method. 
Assuming all non-zero matrix elements have the value 
of oder of 1, the first term in (|3^ ) becomes oder of 
N , while the fiuctuation of the second term becomes 
oder of \/N for sparse matrices since the number of 
non-zero matrix elements is oder of N for sparse ma- 
trices. Therefore, the statistical error of the trace will 
become small as 1/\/N. For example, the statistical 
error becomes 10"'^ for N — 10®, which can be con- 
sidered as small enough. If the statistical error with 
a single set of random variables is not small enough, 
we can repeat the calculation with M sets of random 



variables 



= ,N) where (m = 1, • ■ ■ , M ) 
and obtain the statistical error of order of 

If the operator A is Hermitian, the imaginary part 
of the statistical error becomes zero during the Monte 
Carlo process since 



'>{n'\A\n) 



2Re(e*(^"-^"''(n'|^|n)) 



(35) 



n>n' 



This algorithm can be apphed for evaluating, for 
example, the density of state. 



p[ijj) = — Im G„„(w + i-q) 



— Im (tr [G(a; + i?/)]) . 

TT 



(36) 
(37) 



3.2. Finite Temperature Average of Operators 

The Monte Carlo method for computing trace makes 
it possible to evaluate the expectation value of an ar- 
bitrary operator ^ at a finite temperature T, which is 
defined as 

{A)t - Z-Hr [e-'^ff A] 

= Z-^Y.inle-'^^Aln) (38) 

n 

Z = tr[e-'5«] =^(n|e-''«|n) (39) 

n 

where (3 =1/T (we use fc^ = 1) and the sum runs over 
a chosen basis set of complete orthonormal basis \n). 
In principle, we can compute (|3^) by using the Monte 
Carlo scheme ( |33| ) for evaluating the trace. However, 
the difficulty in evaluating the exponential operator 
g-/3^f n^ay hinder us from applying this straightforward 
scheme. To overcome this obstacle, we transform the 
expression by using the eigenkets of the Hamiltonian 
as a basis set, namely, 

{A)t X Z 

= tr[e-'^^A] =^e-^^^(A|A|A) (40) 

A 

= J2 f dEe-^^{X\S{H -E)A\X) (41) 



X {X\ — iG{E + uj) - G{E - iii))A\X)m 
2i-K 



-1 



dEe 



•KtT:[{G{E + iri)~G{E~iri))A] (43) 

^ _i / dEe-P'^ 
2nrJ 

X (tr [AG{E + ir])] - tr [A'fGiE + ir])] *()44) 
If ^ is a Hermitian operator, ( |4^ ) reduces to 

{A)t X Z 

= — J dEe~^^l-m {ii[AG{E + iri)]) . (45) 

The partition function Z can be evaluated by using the 
unit matrix / in place of A. Note that the imaginary 
part of the Green's functions works as a energy filter 
function extracting the component of energy E from 
the random ket 1$) when we evaluate the trace. 



+ 

e 
e 

^5 



1 






1 




A 


f 


]f 


\ f 


- \ 


















- 












1 


1 1 



-1 
u> 



Fig. 1. ImG„ 



a; - 



irj) for = 10 and rj = 0.1. 



4. Numerical Examples 

In this section, we show several numerical results to 
demonstrate the effectiveness of Particle Source Method. 
For simplicity, we calculate only one-body problems. 
However, these results include the Hamiltonian matri- 
ces of TV = 10^, which is comparable to the dimension 
of the Hamiltonian matrices in manybody problems. 
Therefore we believe that our method is the effective 
also in manybody problems. All numerical results in 
this section have been calculated with complex double 
precision arithmetic of FORTRAN. 

4-.1. Perfect ID Lattice 

Let us study the Hamiltonian of an electron in one 
dimensional space, 



H 



P 

2mp 



(46) 



where rrie is mass of electron and V{x) is the static 
potential. After discretizing in space with the lattice 
size Ax, the Hamiltonian is approximated by a tight 
binding form. 
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n=l 



\ ^ / f t 
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(47) 



where e„ = V{xn) and c]j and c„ are the creation 
and annihilation operator of electron at the site x„ = 
n X Aa; (n = 0, 1, • • • , N). The periodic boundary con- 
dition is set as 

(n = 0|(/)) = {n = iV|(/)) (48) 

where In) is the electron state at the n-th site. 
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Fig. 2. (a) real part and (b) imaginary part of Gnn{'-^ + iv) for ^ ~ 10^ and r] — 10 



When V{x) = 0, the exact analytical eigenstates potential, 
and eigenvalues of the Hamiltonian ( p7| ) with the bound- 
ary condition (Q) are well known, V{x) ~ 
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\E„ 



A exp (ikmnAx) 



E„ 



kn 



n=l 

mn 



[1 — cos (kmAx) 



(50) 

where A is a normalizing constant and m is an inte- 
ger m = 0,±1,±2,---,±(A;'- 2)/2,N/2 for even N 
and m = O, ±1, ±2, • • • , ±(7V - l)/2 for odd iV. Note 
that (^o|) approximates well the parabolic dispersion 
relation ( ^6| ) of the continuum model, if m <C iV or 
E ^ 1. In the following, we set h — nie — Ax — 1 for 
simplicity. 

Figure |l| shows the imaginary part of the Green's 
function G{oj + irj) for = 10 and r] = 0.1, where 
uj = E— 1 is the energy measured from the band center. 
The numerical result reproduces faithfully the exact 
spectrum (5^) of the Hamiltonian (^). 

Figure ^compares the Green's function G„„(aj + 
irj) of a long perfect lattice calculated by the multiple 
frequency method to the exact analytical result. For 
the numerical calculation, we used parameters, a — 
0.1, f] = IQ-^, S = 10^2^ and Acj = 5 x 10"^ and 

= 10^ . The computational time was 3 hours on 
the supercomputer at RIKEN. The exact result in the 
limit N ^ oo and 77 — > -|-0 is calculated by using the 
analytical expression^-*, 



G„„(cj + irj) = 



_sgnM_ 

Vlu2 - 1 



(1^1 < 1) 
(1^1 >1) 



(52) 



4-2. Resonant Scattering by a square well potential 

The transmission probability of a particle described 
by the Hamiltonian (^6|) with an attractive rectangular 
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(49) has an analytical expression 
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Figure ^ compares the transmission probability T^n 
\GlrX(jJ + iri)vf calculated by using PSM with the an- 



alytical result (54). For the numerical calculation, we 
used parameters, a — 0.1, 77 = lO^"*, S = 10^^, and 
Alu = 5 X 10~^ and N — 10^ . The transmission 
probability calculated by PSM is slightly smaller than 
the exact result. This is probably because of the fi- 
nite imaginary part of the energy, rj = 10~*, which 
physically corresponds to the absorbtion of the parti- 
cle. The computational time for this calculation was 3 
hours on the supercomputer at RIKEN. 

5. Summary 

In this article, we developed the PSM for calcu- 
lating the real time Green's functions of large sparse 
N X N Hamiltonian matrices, which exploits the nu- 
merical solution of the inhomogeneous time-dependent 
Schrodinger equation. The method has a clear-cut 
structure reflecting the most naive definition of the 
Green's functions, and is very suitable to parallel and 
vector supercomputers. It requires, as the Lanczos 
method, memory space of oder of N, and the CPU 
time of oder of TV for a given set of parameters. The 
PSM can also calculate matrix elements of the prod- 
ucts of several Green's functions and other operator, 
while the Lanczos method can calculate matrix ele- 
ments of operators that contains only one Green's func- 
tion. This is because PSM can calculate N matrix 
elements, G{u! + ir])\j), at one calculation, while the 
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Fig. 3. Transmission probabil- 

ity Tlh = \Glr{oj + i'ri)vf as a function of energy. 
The parameters are N = 10^ and r] = 10~^. 

Lanczos method can calculate only one diagonal ma- 
trix clement {j\G{uj + i'ri)\j) at a time. We applied the 
PSM to simple lattice models and demonstrated that 
the method can be a powerful tool to study dynamical 
properties of finite quantum systems. 
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